subroutine init_fft()

  use storage
  integer :: i_x,i_y,i_z,is
  
  allocate(k_vec(ns))
  allocate(k1_vec(nt_y/2+1),k2_vec(nt_z))
  
  allocate(spe_vex(ns,nt_x),spe_vey(ns,nt_x),spe_vez(ns,nt_x),spe_vet(ns,nt_x))

  allocate(fou_vex(nt_x,nt_y/2+1,nt_z),fou_vey(nt_x,nt_y/2+1,nt_z),fou_vez(nt_x,nt_y/2+1,nt_z))
  
  do i_y=1,nt_y/2+1
    k1_vec(i_y) =2.*pi/ly*(i_y-1)
    if(i_y .eq. nt_y/2+1) k1_vec(i_y)=2.*pi/ly*(i_y-nt_y-1)
  enddo
  
  do i_z=1,nt_z
     if(i_z .le. nt_z/2) k2_vec(i_z)=2.*pi/lz*(i_z-1)
     if(i_z .gt. nt_z/2) k2_vec(i_z)=2.*pi/lz*(i_z-nt_z-1)
  enddo
  
  k_max = sqrt((2.*pi/ly*(nt_y/2-1))**2+(2.*pi/lz*(nt_z/2))**2)
  dk=k_max/(ns-1)
  do is=1,ns   
     k_vec(is)=(is-1)*dk
  enddo
  
  pos_sci = 0.
  spe_vex = 0.  ;  spe_vey = 0.  ;  spe_vez = 0.  ;  spe_vet = 0.

end subroutine init_fft

subroutine transform_fft()

  use storage
  use iso_c_binding
  include 'fftw3.f03'
!  include 'fftw3.f'

  integer :: i_x,i_y,i_z,is

  do i_x = 1,nt_x
     call sfftw_plan_dft_r2c_2d(plan,nt_y,nt_z,tot_rhu(i_x,1:nt_y,1:nt_z),fou_vex(i_x,1:nt_y/2+1,1:nt_z),FFTW_ESTIMATE)
     call sfftw_execute_dft_r2c(plan, tot_rhu(i_x,1:nt_y,1:nt_z), fou_vex(i_x,1:nt_y/2+1,1:nt_z))
     call sfftw_destroy_plan(plan)
     fou_vex(i_x,1:nt_y/2+1,1:nt_z) = fou_vex(i_x,1:nt_y/2+1,1:nt_z)/(nt_y*nt_z)
   
     call sfftw_plan_dft_r2c_2d(plan,nt_y,nt_z,tot_rhv(i_x,1:nt_y,1:nt_z),fou_vey(i_x,1:nt_y/2+1,1:nt_z),FFTW_ESTIMATE)
     call sfftw_execute_dft_r2c(plan, tot_rhv(i_x,1:nt_y,1:nt_z), fou_vey(i_x,1:nt_y/2+1,1:nt_z))
     call sfftw_destroy_plan(plan)
     fou_vey(i_x,1:nt_y/2+1,1:nt_z) = fou_vey(i_x,1:nt_y/2+1,1:nt_z)/(nt_y*nt_z)
   
     call sfftw_plan_dft_r2c_2d(plan,nt_y,nt_z,tot_rhw(i_x,1:nt_y,1:nt_z),fou_vez(i_x,1:nt_y/2+1,1:nt_z),FFTW_ESTIMATE)
     call sfftw_execute_dft_r2c(plan, tot_rhw(i_x,1:nt_y,1:nt_z), fou_vez(i_x,1:nt_y/2+1,1:nt_z))
     call sfftw_destroy_plan(plan)
     fou_vez(i_x,1:nt_y/2+1,1:nt_z) = fou_vez(i_x,1:nt_y/2+1,1:nt_z)/(nt_y*nt_z)
  enddo

end subroutine transform_fft

subroutine post_fft()

  use storage
  integer :: i_x,i_y,i_z

  do i_x=1,nt_x
     do i_z=1,nt_z
        do i_y=1,nt_y/2+1
           k_mod=sqrt(k1_vec(i_y)**2+k2_vec(i_z)**2)
           if(i_y+i_z .gt. 3) then
             if(i_y .eq. 1) then
                pos_sci(i_x) = pos_sci(i_x) &
                                  + 0.5*(Real(fou_vex(i_x,i_y,i_z))**2+AImag(fou_vex(i_x,i_y,i_z))**2)/k_mod &
                                  + 0.5*(Real(fou_vey(i_x,i_y,i_z))**2+AImag(fou_vey(i_x,i_y,i_z))**2)/k_mod &
                                  + 0.5*(Real(fou_vez(i_x,i_y,i_z))**2+AImag(fou_vez(i_x,i_y,i_z))**2)/k_mod
             else
                pos_sci(i_x) = pos_sci(i_x) &
                                  + (Real(fou_vex(i_x,i_y,i_z))**2+AImag(fou_vex(i_x,i_y,i_z))**2)/k_mod &
                                  + (Real(fou_vey(i_x,i_y,i_z))**2+AImag(fou_vey(i_x,i_y,i_z))**2)/k_mod &
                                  + (Real(fou_vez(i_x,i_y,i_z))**2+AImag(fou_vez(i_x,i_y,i_z))**2)/k_mod
             endif
  !           print*,'vx:',Real(fou_vex(i_x,i_y,i_z)),AImag(fou_vex(i_x,i_y,i_z))
  !           print*,'vy:',Real(fou_vey(i_x,i_y,i_z)),AImag(fou_vey(i_x,i_y,i_z))
  !           print*,'vz:',Real(fou_vez(i_x,i_y,i_z)),AImag(fou_vez(i_x,i_y,i_z))
  !           print*,'k_mod:',k_mod
  !           print*,'premi tasto please' ; read(*,*)
           endif
        enddo
     enddo
  enddo

  do i_x=1,nt_x
     do is=1,ns
        do i_z=1,nt_z
           do i_y=1,nt_y/2+1
              k_mod=sqrt(k1_vec(i_y)**2+k2_vec(i_z)**2)
              if((k_mod.lt.(k_vec(is)+dk/2.)).and.(k_mod.gt.(k_vec(is)-dk/2.)))then
                if(i_y .eq. 1) then
                  spe_vex(is,i_x) = spe_vex(is,i_x) + 0.5*(Real(fou_vex(i_x,i_y,i_z))**2+AImag(fou_vex(i_x,i_y,i_z))**2)/dk
                  spe_vey(is,i_x) = spe_vey(is,i_x) + 0.5*(Real(fou_vey(i_x,i_y,i_z))**2+AImag(fou_vey(i_x,i_y,i_z))**2)/dk
                  spe_vez(is,i_x) = spe_vez(is,i_x) + 0.5*(Real(fou_vez(i_x,i_y,i_z))**2+AImag(fou_vez(i_x,i_y,i_z))**2)/dk
                else
                  spe_vex(is,i_x) = spe_vex(is,i_x) + (Real(fou_vex(i_x,i_y,i_z))**2+AImag(fou_vex(i_x,i_y,i_z))**2)/dk
                  spe_vey(is,i_x) = spe_vey(is,i_x) + (Real(fou_vey(i_x,i_y,i_z))**2+AImag(fou_vey(i_x,i_y,i_z))**2)/dk
                  spe_vez(is,i_x) = spe_vez(is,i_x) + (Real(fou_vez(i_x,i_y,i_z))**2+AImag(fou_vez(i_x,i_y,i_z))**2)/dk
                endif
              endif
           enddo
        enddo
     enddo 
  enddo

  do is=1,ns
     spe_vet(is,1:nt_x) = (spe_vex(is,1:nt_x) + spe_vey(is,1:nt_x) + spe_vez(is,1:nt_x)) 
!     print*,'is,spe_vet(is):  ',is,spe_vet(is)
  enddo

end subroutine post_fft

subroutine reduce_fft()

  use storage 

  pos_sci(1:nt_x) = pos_sci(1:nt_x) / ((post_itmax-post_itmin)/post_itout+1)
  do is=1,ns
     spe_vet(is,1:nt_x) = spe_vet(is,1:nt_x) / ((post_itmax-post_itmin)/post_itout+1)
!     print*,'is,spe_vet(is):  ',is,spe_vet(is)
  enddo

end subroutine reduce_fft

subroutine print_fft()

  use storage 
  integer :: is
  real(myk), parameter :: ni=0.00005_myk
  real(myk), dimension(1:nt_x) :: l,u_rms,Re_l,l_k,l_o,Re_o,t_o

  open(unit=10,file='spettro.dat')
  print*,'ns: ',ns
  do i_x=1,nt_x
     do is=1,ns
        write(10,*) ix,is,spe_vet(is,i_x)
     enddo
  enddo
  close(10)

  u_rms = sqrt(2.*pos_ene/3)
  l_o   = (pi/2.)*pos_sci/u_rms**2
  t_o   = l_o/u_rms
  l     = sqrt(5.*pos_ene/pos_ens)         
  Re_l  = u_rms*l/ni
  l_k   = (ni**2/(2.*pos_ens))**.25
  Re_o  = u_rms*l_o/ni

  open(unit=10,file='parametri.dat')
     write(10,*) 'u_rms, Scala integrale,Scala di Taylor,Reynolds lambda,Reynolds zero,Eddy turn-over time,Scala di Klomogorov '
     do i_x=1,nt_x
        write(10,'(100(F20.12,1X))') u_rms(i_x), l_o(i_x), l(i_x), Re_l(i_x), Re_o(i_x), t_o(i_x), l_k(i_x)
     enddo
  close(10)

end subroutine print_fft
